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Abstract 

> 

0^ I A new classical interaction potential for water simulations is presented. Water 

is modeled as a fully dissociable set of atoms with a point dipole, determined self- 
, consistently, on every oxygen atom. The oxygen polarizability is not fixed but de- 

I pends on the geometry of the system. We show that, in spite of the limited number 

of free parameters, the model reproduces the geometrical and vibrational properties 
of microclusters in a satisfactory way. 



S 

^ I Water, being arguably the most important liquid in our life, has always at- 

O ■ tracted a great deal of attention from physicists and chemists. It is also a rather 

difficult substance to model, and, in spite of vast improvements in numerical 
techniques, ab-initio simulations are not yet able to simulate more than a 
^ I hundred molecules for a reasonable time. Many classical potentials have been 

cd proposed, but very few of them are able to take into account all the important 

structural and dynamical features of water. Dissociability, for instance, is im- 
portant in order to describe bulk water and its hydration and proton-transfer 
properties. Polarizability is crucial in the description of structural properties. 

In this Letter we propose a new dissociable and polarizable potential for wa- 
ter and validate it by means of the study of the geometrical and vibrational 
properties of the water molecule and of the water dimer and trimer; results on 
larger clusters, and a full comparison with other potentials, will be presented 
in a full paper. The validation process we used is especially stringent as we 
predict a range of physical properties much larger than usually done in this 
kind of studies. We believe this gives a better guarantee of transferability to 
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our potential than the usual procedures. Our potential contains a small num- 
ber of parameters and is based on reasonable physical assumptions, which 
include some quantum effects previously neglected. Its aim is to reunite the 
advantages of dissociable and polarizable potentials while at the same time 
keeping a comparatively simple and physically transparent form. 

We assume that, in the molecule, most of the electric charge of each hydrogen 
is transferred to the oxygen. We also assume that the electronic polarizability 
of the molecule is entirely due to the oxygen. Several models already exist that 
allow one to treat successfully ionic crystal structures (see for instance [1,2,3]), 
and we will start from there. The potential we propose has the following form: 

U = Uzz + Ubm + Udd + Uvw + UzD + Usr. 

In the term 

C.. = iE^ (1) 



Zo and Zh indicate ionic charges and are considered as adjustable parameter, 
bound by the charge neutrality constraint Zo -\- 2Zh — and Vij — Ti — rj, 
where Tj is the position of the i-th nucleus. The Born-Mayer repulsive term 

Ubm-W f{p. + P,)exp f ^^^^^^) (2) 



includes the adjustable parameters /, i?o, Rhi Poi Ph- The term 

Udd = -2_^ 2^ —3 -5 ]+2^7^ i^) 




describes dipole-dipole interaction, Pj being a point dipole located on the i— th 
oxygen, and the bare oxygen polarizability cvo is taken as a further parame- 
ter; since hydrogen is almost entirely deprived of electrons, it is considered 
unpolarizable. The Van der Waals term is written as 

Uvw = -\Y.Y."^ (4) 

^ i j¥=i 'ii 



while the term 

Uzn = Y.T. P. • Z-f (5) 
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describes the charge-dipole interaction. The term 

ieH jeo '1-3 '-"■o Pi- Pj 



is taken from Ref.[2], where the polarizabihty in the shell model is studied. In- 
tuitively, this term reduces the polarizabihty of oxygen by taking into account 
the interaction between its electronic charge and that of the closest hydrogen 
atoms, so it is not surprising that it has the same form of Eq. 2. In principle 
the values considered for Rsr„ + Rsm (o^ily the sum of these parameters is 
relevant) are independent of Rg + Rh, but in practice they are only slightly 
different, and in some previous studies [2] they were set as equal. 

In almost all this work we set ph = 0, that is we neglected the short-range 
hydrogen-hydrogen interactions given by Eqs. 2 and 6. In short, our potential 
is characterized by 11 independent parameters, which are listed in table 1 

together with their o])timized values. 



Zh = 0.80e 


Ph 


= 


Zo = -2Zh = -1.60e 


Po 


= 0.2091 


/ = 0.03405e2l-2 




= 1.0271''^ 


Rh = 0.0121 


F 


= 0.296-^1^ 


Ro = 1.7181 


Ch 


= O.Olel^/^ 




Co 


= 3.9el5/2 



Table 1 
Optimized set of parameters 

It is worth noting that the optimized value of cto lies approximately midway 
between the polarizabilities of oxygen and of the water molecule. We fitted 
several physical quantities we considered of primary interest: of the monomer, 
the geometrical data, the vibrational frequencies and the dipole moment, all 
of them known experimentally. Of the clusters, the geometrical data relative 
to the equilibrium configurations, which are known experimentally and some- 
times from ab-initio calculations. The equilibrium geometries were found us- 
ing a simple steepest-descent algorithm, with determined self-consistently, 
while the vibrational frequencies were obtained via numerical calculation and 
diagonalization of the dynamical matrix. After fixing the parameters of the 
potential, as a further test we calculated the cluster binding energies, the vi- 
brational frequencies and some configurations considered as local minima of 
the potential energy as reported also in some ab-initio studies, comparing all 
these data with those found in the literature. 

In the rest of this paper we shall use the following symbols: ii for the dipole mo- 
ment of the molecule; Vss, ^as and z/'b for the molecule normal modes, that is, re- 
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spectively, symmetric stretching, asymmetric stretching and bending; Jlimolec) 
is the average molecular dipole moment in the clusters. The subscripts h and 
/ indicate the bond and free hydrogens, respectively. 





Exptl data 


ab-imtio[14] 


This work 


d(O-H) 


0.957A[4] 


0.965A 


0.961A 


Oh-o-h 


104.5° [4] 


103.8° 


104.5° 


i^h 


1595cm-i[5] 


1623cm-i 


1158cm"^ 




3657cm-i[5] 


3807cm" 1 


3609cm"i 




3756cm-i[5] 


3936cm" ^ 


3234cm"i 


^J' 


1.86D[6] 


1.87D[9] 


2.07D 



Table 2 
Molecule data 



The results for the molecule are shown in table 2. As one can sec, the overall 
agreement with the experimental data is good, the only quantity in error of 
more than 15% being v^. The frequencies i/ss and i/as are in the wrong order, a 
feature common in our potential, but given the small difference between them 
this does not look too worrying. Attempts made to fix the discrepancy led to 
an oxygen-oxygen repulsion that is too small. 

The numerical results on the dimer are collected in table 3. The ground-state 
geometry is, qualitatively, the usual one; is the angle between the bisector 
of the H-O-H angle of the donor molecule and the 0-0 direction, 6a is the 
same for the acceptor. is the dimer dissociation energy. 





Exptl* / ab-initio^ data 


This work 


d{0 - O) 


2.98A[7]* 


3.20A 


B{0 -H^-O) 


174° ±10° [7]* 


171° 


d{0 - H)b 


0.972A[14]t 


0.98A 


d{0-H)f 


0.964-0.966A[14]t 


0.960-0.968A 


0a 


57° ± 10° [8]* 


92° 




51° ± 10°[8]* 


57° 


De 


0.23eV[7]* 


0.16eV 


Jl{molec) 


2.10D[9]t 


2.22D 



Table 3 

Dimer structural data 



We found that, as Co and Ro vary, when the oxygen-oxygen separation increases 
then ao decrease by about the same amount, and vice versa: this is clearly 
an electrostatic effect due to the attraction of the hydrogens in the acceptor 
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molecule and the oxygen of the donor molecule. In addition, the angle 9a 
increases, as one would expect, as the molecular dipole moments decrease and 
tend not to align. The excessive value of 6a is therefore a price we have to 
pay in order to obtain a dipole moment not too high. We also verified that, as 
d{0 — 0) increases, decreases. It is thus not surprising that, since d{0 — 0) 
is overestimated, De is underestimated. 





H/XpZl Uillil 


Ab-initio (MP2) 


1 nib W07 rv 


Intermolecular 


[10] 


[14] 


[11] 




LV12: Donor torsion 


- 


141 


155 


107 


LOii: Acceptor twist 


- 


147 


193 


238 


a;8: Acceptor bend 


243 


155 


178 


262 


UJ7: 0-0 stretch 


155 


185 


220 


94 


uq: In-plane donor wag 


320 


342 


398 


487 


uJiQ-. Out of plane donor bend 


520 


632 


715 


722 


Intramolecular 


[12] 








H-O-H acceptor bend 


1601 


1624 


1615 


1196 


H-O-H donor bend 


1619 


1651 


1636 


1303 


Donor bridge 0-H stretch 


3530 


3712 


3539 


3127 


Acceptor sym. 0-H stretch 


3600 


3799 


3659 


3590 


Donor free 0-H stretch 


3730 


3907 


3719 


3550 


Acceptor asym. 0-H stretch 


3745 


3929 


3746 


3215 



Table 4 



Dimer normal modes (in cm ^) 

Result for the dimer normal modes are shown in table 4. For their definition 
we follow Ref.[14] and their identification Refs [10,20,11]. The ordering of the 
intermolecular frequencies we find is not in complete agreement with experi- 
mental or ab initio data; it is the same found by [14], with the exception of 
the "0-0 stretch" mode. The experimental data show that the intramolecular 
frequencies are not much different from those of the monomer; the main dif- 
ference consists of an increase of the bending frequencies and a decrease of the 
symmetric and asymmetric stretching frequencies. These trends, confirmed by 
ab initio calculations ([14] and [19]), are well reproduced in our results. 

The equilibrium geometry of the trimer is shown in fig.l and the numerical 
data are collected in table 5. The longest side of the triangle corresponds to the 
smallest H-bond angle. This is characteristic of our potential. We obtained two 
relative minima, shown in fig. 2 and 3. The first is 0.06 eV above the energy 
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Fig. 1. Trimer equilibrium geometry 




Fig. 2. Trimer: the lowest metastable state 

in the absolute minimum, the second 0.15 eV. Their ordering in energy is the 
same as in the ab initio calculations [17,18], respectively 0.04 and 0.24 eV. 

As an illustration of the overall quality of our potential, in Table 6 we show 
the trends of two important quantities — the average oxygen-oxygen distance 
and the average molecular dipole moment — with increasing number n of 
molecules in the cluster. To make trends clearer, for n < 8 the lowest-energy 
cyclic clusters have been considered, even when they are not the absolute 
energy minima. It is clear that the trends from ab initio calculations are well 
reproduced and the absolute errors actually decrease with n. 

In conclusion, we have shown that the model reproduces the basic physical 
properties of microclusters in a satisfactory way. In particular, we are able 
to obtain bond lengths and dipole moments within 10-20% of experimental 
and ab-initio values, and to reproduce qualitatively not only all the ground- 
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ijxp. acLvd 


1 Ihlb WOftv 


GL[KJ-KJ ) 


O/l Q7 07 A 


0. iu-o.uy-o.zoA 


U{Ui — tlh — U2) 




IDo 


a(y2 — Hh — Us) 


iOU 


iOD 


6/(C3 - Hb-Ui) 


ICO 


164 


TT ID /n \ 
tih\ - (Ci — C2 — t^a j 


ioU 


i (0 


-"62 — {y2 — C3 — Cl j 


1 sjn 
ioU 


1 s2n 

ioU 


J^b3 i*^3 Ui — U2) 


1 sn 

IOU 


1 81 
1(51 




UjU itllilU 




ri(n n\, 
a[u — ri )i) 


n Q77A<^ 

u.y / 1 -TV 


n QQ 1 nnA 
u.yy-i.uuri. 


OiXkJ — ) f 




u.yDo-u.yoo-u.yoD-tA. 


-n/i — \Ui — C2 — Lv'sj 


iio 


yo 


-H/2 (,t>'2 — — Ui) 


9076 


zoo 


-n/3 — i<^3 — (^2 — ^1) 


Oqi e 


zoo 




U.zzeV 


U.loeV 


Jiirnolec) 


Z.61U 


2.61 u 


Frequency range (cm^^) 






intermolecular 


1 ni Tori b 


/0-lUoo 




ICO or'O G 
LOO-006 






1 Qc nci f 




intramolecular 


1753-4215 ^ 


1288-3544 




1632-3898 " 






1623-3727 ^ 





Table 5 

Trimer data "Ref.[15].^Ref.[14], HF.'=Ref.[16], LDA+GC.''Ref.[9].'=Ref.[14], 
MP2/Ref.[19], LDA+GC. 

state geometries, but also several local minima. This remains true when the 
potential is applied to larger clusters, as we will show in a forthcoming paper. It 
should be stressed that for almost none of the potentials found in the hterature 
one can find a comparison with experimental or ab-initio data as extensive 
and thorough as ours. Moreover, the limited number of parameters we used 
makes it more likely that their values are close to the real ones. We therefore 
believe that our potential is a promising one with regard to transferability and 
accuracy. 
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Fig. 3. Trimer: the next-lowest metastable state 



n 


1 


2 


3 


4 


6 


8 


d{0 - 0){A) This work 




3.20 


3.14 


2.80 


2.77 




d{0-0){A) [14] 




3.03 


2.93 


2.88 


2.86 




Jl{molec){D) This work 


2.07 


2.22 


2.37 


2.76 


2.82 


2.77 


Ji{molec){D) [9] 


1.86 


2.10 


2.31 


2.55 


2.70 


2.72 



Table 6 

Trends with varying n in cyclic clusters and the octamer, compared with ah initio 
results 

This work has been funded by a MIUR-COFIN grant. 
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